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Abstract 

We discuss various critical densities in sandpile models. The sta- 
^vq tionary density is the average expected height in the stationary state 

of a finite-volume model; the transition density is the critical point 
(-h in the infinite- volume counterpart. These two critical densities were 

Qh generally assumed to be equal, but this has turned out to be wrong 

<— | for deterministic sandpile models. We show they are not equal in a 

quenched version of the Manna sandpile model either. 

In the literature, when the transition density is simulated, it is 
i— i implicitly or explicitly assumed to be equal to either the so-called 

| threshold density or the so-called critical activity density. We properly 

J> define these auxiliary densities, and prove that in certain cases, the 

threshold density is equal to the transition density. We extend the 
^ definition of the critical activity density to infinite volume, and prove 

■^j- that in the standard infinite volume sandpile, it is equal to 1. Our 

results should bring some order in the precise relations between the 
various densities. 



Keywords: Sandpile model, phase transition, critical density, transition 
density, Manna model, activity density. 



^ AMS Subject Classification: 60K35, 60K37, 60J05, 60G10 



1 Introduction 

This paper aims to bring some order in a long-standing discussion on the 
equality or inequality of various critical densities in sandpile models. By 
means of introduction and context, we first informally describe the archety- 
pal Abelian sandpile model (ASM or 'classical' sandpile model). On I n = 
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{0, 1, . . . , n} we initially assign a height (a non- negative integer) to each site 
x G /„. The dynamics is as follows. Each time step starts with an addition: 
the height of a randomly chosen site is increased by 1. Then the configu- 
ration is stabilized by so called topplings: a site with height 2 or more is 
unstable and may topple. In a toppling, its height decreases by 2 and the 
height of both neighbors increases by 1. Sites or n have only one neighbor; 
when these sites topple, one particle leaves the system. The configuration is 
stabilized when there are no more unstable sites to topple. It is not hard to 
show that in the stationary state of this model, the expected height is 1 — 
for every site; see for example [24, 29J. We call this expected height, in the 
limit as n — > oo, the stationary density and denote it by p s . (So in this case, 

Ps = I-) 

Next consider a related model, the infinite volume model. On the infinite 
line Z, start with (independent) Poisson-p number of particles at each site. 
At each discrete time step, we simply topple all unstable sites simultane- 
ously. (There are alternatives, like toppling after exponential waiting times, 
independently for each site, see [19] for a detailed discussion of this.) It is 
known [TH] that for p < 1, only finitely many topplings per site are needed 
to obtain at a limiting stable configuration, whereas for p > 1, every site 
topples infinitely many times. Hence, this model has a phase transition at 
p = 1 and we call this the transition density p tr - 

In a widely cited series of papers [H |2|, [3J, HI |5] , Dickman, Munoz, Vespig- 
nani and Zapperi developed a theory of self-organized criticality (SOC) as a 
relationship between driven dissipative systems like the ASM, and systems 
with particle conservation like the infinite volume model. Loosely speaking, 
the following heuristic connection (the 'density argument') is offered: in the 
presence of activity of the ASM, that is, when topplings take place, we should 
typically have p > p tr , and particles are leaving the system. In the absence of 
activity, we typically have p < p tr , and there is addition of particles. Hence 
the density of particles always changes in the direction of p tr and therefore 
the only candidate- value for the stationary density is pt r - 

This would be a very attractive explanation for the observation that the 
stationary state of the two-dimensional ASM shows power law behaviour of 
many spatial and temporal quantities, (apparently) without any tuning of 
a parameter. However, it requires that p s = p tr . In the one-dimensional 
case this is true. At first, simulations and exact calculations for the two- 
dimensional ASM seemed to confirm that p s and pt r are equal up to three 
decimals. However, more accurate simulations show a difference in the fourth 
decimal for the two densities on Z 2 . Moreover, on several other graphs they 
are provably different [TS]. It is at present believed that this difference is a 
consequence of the existence of so called toppling invariants [UJ • 
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The more recent literature (we give references in Section [4]) on the density 
argument tends to focus on the Manna model, a sandpile model with stochas- 
tic topplings, for which it is believed that toppling invariants play no role. 
The Manna model is defined as above, with the important difference that 
during a toppling, particles go to randomly chosen neighbors, independently 
of each other, and independent of earlier topplings of the same or other sites. 
In other words, the topplings in the Manna model are subject to annealed 
disorder. 

Simulating p s for a sandpile model is straightforward, since this density is 
defined for a model in finite volume. However, p tT is defined in infinite volume, 
which makes simulations much more difficult. In reality, when people try to 
simulate p tr , what they really simulate is the so called critical activity density 
or the threshold density, which we define and study below. It is then tacitly 
assumed (sometimes also explicitly) that these densities are equal to p tr - 

The goal of the present paper is to carefully distinguish between the var- 
ious densities, and to investigate the density conjecture in quenched versions 
of the Manna model, that is, in sandpile models in which the toppling rules 
per site are random but fixed throughout the evolution of the process. It 
turns out that we can explicitly compute p s and p tr in these quenched Manna 
models, and that in general they are not equal. We also study the relation 
between the various densities in the original (annealed) Manna model. 

We end this introduction with a precise definition of the sandpile models 
which we consider. We distinguish between 

1. Driven models, that is, models with addition of particles and topplings, 
always on the interval = {0,1, ... ,n} with dissipation at the open 
boundaries; 

2. Infinite volume models on Z, without additions, and only topplings; 

3. Fixed-energy models, that is, models with topplings but without addi- 
tions, on the n-cycle C n — {0, 2, . . . , n}, where and n are neighbors. 

We abuse notation and view I n , C n and Z as graphs with nearest neighbor 
connections. A typical configuration will be denoted rj and for any vertex 
(site) x, r)(x) G {0, 1, 2, . . .} is interpreted as the height at site x. We call a 
site stable if its height is or 1. If the height of a site is 2 or more, then we call 
this site unstable. If all sites are stable, then we say that the configuration 
is stable. 

Upon an addition at site x, the height of site x increases by 1, and can 
henceforth become unstable. An unstable site may topple. With a toppling 
of site x, the height of site x decreases by 2, and the height of the neighboring 
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sites may increase, in such a way that the total increase is at most 2. We 
denote the neighbors of site x as x — 1 and x + 1, where on the graph C n , 
we interpret this modulo n + 1. On J n , and n only have 1 neighbor. 
We consider three possibilities for the height increase of the neighbors in a 
toppling, namely 

1. The height of site x — 1 (if present) increases by 2, and the height of 
site x + 1 does not change. We call this a left toppling (I). 

2. The height of site x + 1 (if present) increases by 2, and the height of 
site x — 1 does not change. We call this a right toppling (r). 

3. The height of both sites x — 1 and x + 1 (if present) increase by 1. We 
call this a symmetric toppling (s). 

In this paper, we introduce quenched sandpile models, which we define 
as follows. On a graph G, we let v be a probability measure on {/,r, s} G , 
an element of which we call a toppling environment. We choose a toppling 
environment £ according to v. Given a toppling environment (, we allow 
at site x only topplings of type ((x). We say that ((x) is the Za6eZ at site 
x. The combination of a configuration 77 and a toppling environment ( is 
called a siaie, and denoted (77, £)• We call this state stable if 77 is a stable 
configuration. 

If a state is unstable then topplings can lead to new unstable sites. We 
say we stabilize a state (77, £) if we keep toppling unstable sites, according to 
the (fixed) toppling environment, until there are no unstable sites anymore. 
It may happen though that we never obtain a stable configuration this way; 
in this case we say that the state is not stabilizable. If a stable configuration 
can be reached, then the order in which we perform topplings has no effect 
on the final stable configuration; this is the famous Abelian property of the 
sandpile (see [TH] for some technical details on this fact in the case of infinite 
volume). When I- and r- labels are not allowed, then this model reduces to 
the classical sandpile model. 

The Manna model is the annealed version of the above quenched sandpile 
model. One convenient way to describe the Manna model [30J , is that instead 
of a single label, we have a stack of labels at every site. Every stack consists 
of an infinite sequence of i.i.d. random labels, which are s with probability 
1/2, I with probability 1/4 and r with probability 1/4, independent of all 
other stacks. When a site topples, it topples according to the first label in 
its stack, and then removes that label from the stack. 

This paper is organized as follows: In Section [2j we define and state our 
results on the two main critical densities, namely the stationary and the 
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transition density. In Section [3j we discuss the additional critical densities 
that have been introduced in the literature to simulate p tr . We state a number 
of new results on these densities. Finally, we give the proofs of all our results 
in a separate Section [5} 

2 p s and p tr for the quenched sandpile models 

2.1 Driven sandpiles on I n and the stationary density 

Driven sandpile models are defined on I n . Let p be a probability measure on 
{0, 1} N and v a translation-invariant probability measure on {/,r, s} N . We 
use p to draw an initial configuration, and v to draw a toppling environment. 
More precisely, to implement the process on I n , the initial configuration is 
the restriction of a drawing from p to the first n + 1 coordinates, and likewise, 
the toppling environment ( n is the restriction to the first n + 1 coordinates 
of a drawing from v. 

Dynamics is as follows. Given n°, we add a particle at a randomly cho- 
sen site (uniformly over all sites) and stabilize the resulting configuration. 
The result is denoted n\. Repeating this process leads to to a sequence of 
configurations (^), t = 0, 1, 

We are interested in the 'average' number of occupied sites in the process 
on /„ in the limit as n — > oo. For given n and (, the process is a Markov 
chain, but it turns out that for certain toppling environments, this Markov 
chain is periodic so that we do not have a well defined unique stationary 
distribution. Nevertheless, the 'average' hinted at above often exists in an 
unambiguous way. We denote by p l n = v) the expected fraction of sites 
in I n whose height is 1 at time t. Here, expectation is with respect to p x v. 
We say that v is directed if either ^(C(0) = r) = or f (C(0) = I) = 0. 

Theorem 1. Suppose that v is directed. Let f r = ^(C(0) ^ s). Then 

Ps = Ps(p, v) := lim lim p l n (2.1) 

n— »oo t— >oo 

exists and is equal to 




In particular, p s does not depend on p. Note that in the classical case, f r = 0, 
and Theorem Ul reduces to the well known result which we described in the 
introduction. 

The quantity p s is called the stationary density. 
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2.2 Infinite volume sandpiles on Z and the transition 
density 

Infinite volume sandpiles in this article are defined on Z. We choose the 
initial state r\ according to /x p , a Poisson product measure with parameter 
p. The toppling environment ( is chosen as before, according to a stationary 
and ergodic probability measure v on {I, s, r} z . There are no additions in 
this model, only topplings. At any discrete time t, we simply topple all 
unstable sites once. This implies that for every t, the distribution of rf is 
translation-invariant. 

We define the transition density as 

Ptr = Ptriy) = sup{p : (77, C) is stabilizable (fi p x v) — a.s.}. 

As mentioned in the introduction, it was recently proved for the classical 
sandpile model on the bracelet graph [18] (with proper interpretations of p s 
and p tr ) that p s 7^ p tr . The current belief is that this inequality holds in 
general for classical sandpile models, and that the one-dimensional case is an 
exception. 

In this paper, we prove that p s 7^ p tr for quenched sandpiles, as long as v 
is directed and satisfies 

u(( : (( x ) = CO + 1) = s) = and v(((x) = ((x + 2) = s) = 0. (2.2) 

for all x. In words, this condition means that between every pair of s labels 
there must be at least two r labels. Note that the transition density is 
not interesting for non-directed models: indeed, between r- and an Z-labelled 
sites, particles cannot escape, and since it is possible with positive probability 
that too many particles are trapped in between, the model will not stabilize 
a.s., for any value of p > 0. 



Theorem 2. Suppose that v is directed and satisfies (2.2). Let f r = ^(C(0) 
r). Then we have 



(a) p s = l- 



(b) p tT is the unique solution of p — \(1 + e 2p ) = ^-^(l + e 2p ) 2 . 



In particular, p s 7^ p tr for v satisfying (2.2) 
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3 Simulating p tr : more critical densities 



3.1 Fixed-energy sandpile models and the threshold 
density 

The fixed-energy sandpile models in this article are defined on C n . They 
play an important role in simulating p tr . In the initial state, r] n is chosen 
according to p p , a Poisson product measure with parameter p, and ( n is 
chosen according to v. The fixed-energy model evolves in time by topplings. 
If the initial state is stabilizable, then every site topples only finitely many 
times. 

We define the threshold density as 

pth = Pth{v) = lim sup{p : (p p x u)((r] n ,(n) is stabilizable) > 1/2}. 

For the classical fixed-energy sandpile model, the threshold density is equal 
to 1. 

Note that the definition of pth is very similar to that of pt r . It is often 
tacitly assumed that they are equal for the classical sandpile on general 
graphs: in [IS], this assumption is explicitly stated. In the literature, only 
for the classical sandpile in dimension 1 a proof is available. Note that for 
a quenched sandpile model, if both /-sites and r-sites exist, pth = 0, since it 
is always possible to have too many particles between a r-site and an /-site, 
from which particles can never escape. Hence, again, only directed models 
are of interest. We prove the following. 



Theorem 3. Suppose that v is directed and satisfies (2.2). Then p tr = p th 



3.2 Parallel chip firing and the activity density 

Consider fixed-energy models with parallel toppling order. These cellular 
automata are called parallel chip firing models [2SJ. We are interested in the 
activity density p a (p p , v, n), which loosely formulated, is the eventual average 
fraction of unstable sites when parallel chip firing is run on a finite graph, 
C n in our case. The next step is to take the limit as n — > oo: 

p a (Pp,v) = limsupp a (ju p , u,n). 

n— >oo 

For the classical sandpile on C n , the graph of p a versus p has been shown 
[B] to take the shape of a 'staircase' with three stairs (see Figure [T]). More 
specifically, the classical model on C n behaves in the following particularly 
simple way. If the total number of particles is strictly in (n, 2n), then every 
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Figure 1: The staircase graph for the classical sandpile model on C n . The 
horizontal axis is p, the vertical axis is the activity density. The graph is the 
same for the fixed-energy and the infinite volume model; see [0] and Theorem 
® 

initial configuration will eventually have period 2. Thus, as t — > oo, the total 
number of unstable sites in one time step will eventually be close to n/2. 
If the total number of particles is larger than 2n, then eventually all sites 
will be unstable. Finally, it was already known that if the total number of 
particles is smaller than n, then every initial configuration stabilizes, so that 
eventually the total number of unstable sites is 0. Thus, the staircase of the 
classical sandpile model on C n has three stairs, and it does not even depend 
on the precise initial distribution of the particles. Only if the total number 
of particles is equal to n or 2n, there are multiple possibilities for the limiting 
fraction of unstable sites. 

For the flower graph [18], the corresponding plot has five stairs. In the 
case of the two-dimensional torus, it is conjectured, based on numerical obser- 
vations [8], that the plot is a devil's staircase. Levine proved the occurrence 
of a devil's staircase for parallel chip firing on the complete graph, with a spe- 
cial class of initial configurations [2T] . To show that the plot need not always 
be a staircase, we include a simulation for the sandpile model of Theorem [2] 
(see Figure [2]). 

The activity density has been used to simulate p tr , under the assumption 
that p tr is equal to the so called critical activity density p ca defined as 

Pea = Pca{v) ■= SUp{p I p a (p P , v) = 0}. 

One can prove this for a number of specific graphs that have an infinite 
volume counterpart, for instance C n (where both densities are trivially equal 
to 1), or, less obvious, for the bracelet graph [IS]. The flower graph or the 
complete graph however do not have an infinite volume counterpart. 

There is an infinite-volume analog of the activity density. Let a p (n,t) be 



8 



0.8 




Figure 2: Simulation of the the quenched restricted model on iioocb with 
f a = 1/4. There are 10 trials for each of 50 values of p. In each trial, a new 
random initial state was generated. Numerically, the predicted value for p tr 
(Theorem [2]) in this case is 0.556... 

the fraction of active sites in [— n, n] after t iterations of the process, starting 
with p p . For fixed t, the limit a p {t) := lim^oo a p (n, t) exists by ergodicity, 
and s. constant. The following result extends [B] to infinite volume. 

Theorem 4. In the infinite volume classical sandpile model on Z we have 
the following. 

1. For p < 1, lim^oo a p (t) = 0. 

2. For p > 2, lim^oo a p (t) = 1. 

3. For 1 < p < 2, Hindoo i ^f =1 a p {t) = |. 

In words, below p = 1, the fraction of unstable sites converges to 0, above 
p = 2 it converges to 1, and in between the fraction converges, in a weak 
sense, to 1/2. The infinite- volume analog of the critical activity density is 
therefore equal to 1, and in the classical sandpile in one dimension, all critical 
densities are equal. 

4 Critical densities for the Manna model 

It is clear how the stationary density should be defined for the Manna model, 
and it is conjectured that it is strictly less than 1: Sadhu et al [21] claim 
that for the Manna model, p s = 0.953 . . ., based on extrapolation of exact 
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calculations for n < 12. Theorem [T] perhaps supports the plausibility of the 
conjecture that p s < 1 for the Manna model. 

As far clS Pi r IS concerned, Dickman et al. [T7j have estimated p tr for the 
Manna model as 0.94887, in fact estimating the critical activity density on C n . 
Huynh et al. [20] performed extensive simulations for the Manna model on 
many one- and two-dimensional lattices. Their estimates confirm the result 
of [IT], and in general give a value close to 0.9 for one- dimensional lattices, 
and close to 0.7 for two-dimensional lattices (note that the table with particle 
densities is not included in the first version of this paper). Sidoravicius and 
Rolla [30] proved that p tr > 1/4. Sadhu et al. [31J compare their value for p s 
of the Manna model to the estimate for p tr of [17] , implying that they should 
be equal (See also Conjecture [6] below). 

It is not hard to show that for the Manna model, p t h = 1- For complete- 
ness, we state and prove this now. 

Proposition 5. For the Manna model on 7L, p th = 1. 

Proof. We make use of the following fact for the classical sandpile model on 
C n : If the total number of particles is less than n, then the model will stabilize 
after finitely many topplings per site, irrespective of the initial position of 
the particles. This can be deduced in the following manner: if the model 
stabilizes, then at least one site does not topple. Then the neighbors of that 
site topple at most once, and so on. 

Now for the Manna model on C n , we choose p < 1 — e, with e > 
arbitrarily small. Then the probability that the total number of particles is 
more than n(l — e/2), tends to as n — > oo. If the total number of particles is 
less than n, then the model stabilizes a.s. To see this, note that one possibility 
for stabilizing is if there is some time t such that starting from t, all unstable 
sites topple according to an s label until the model stabilizes. Since all labels 
are i.i.d. and only finitely many s labels are needed, this event occurs a.s. 
Therefore, the Manna model on C n stabilizes with probability tending to 1 
as n — 7- oo for all p < 1, so that p t ^ = 1. □ 

We conclude that the critical activity density should be at least 1. How- 
ever, the simulations for the Manna model arrive at a value strictly less than 
1. The discrepancy results from a subtlety in choosing the observation time 
window. It is observed that there is a density value strictly smaller than 1, 
such that for p below this value, the total number of topplings before stabi- 
lization is relatively small, whereas for p above this value, the total number 
of topplings is much larger. It is this value that is taken to be an estimate 



for p tr , rather than p th as defined in (3.1). In the second case, before even- 



tually stabilizing, the models is observed to settle into a 'quasi-stationary' 
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state, where the fraction of active sites is approximately constant for a long 
time; see [16], or Figure 3 of [13] for numerical simulations of this state. It 
is believed that this quasi-stationary state of the fixed-energy Manna model 
reflects the true stationary state of the infinite volume Manna model. There- 
fore, the observation time window is chosen such that the model is observed 
in the quasi-stationary state, rather than large enough so that the model has 
stabilized. 

Thus the findings of [T7] and [H] indicate the following conjecture for the 
Manna model, which we now state explicitly. We stress that the conjectured 
relation between the three densities is markedly different from what is known 
for other models. For instance, for the classical sandpile model on the bracelet 
graph it is known [18] that p s ^ p tr = pth, and combining Theorems [2] and 
[3] of this paper gives the same relation for the directed quenched sandpile 



models satisfying (2.2). 



Conjecture 6. For the Manna model on 7L, 

Ps = Ptr 7^ Pth{= !)■ 



With kind permission of Lionel Levine, we give an example to illustrate 
quasi-stationarity in the Manna model. In this model, there is a value of p 
such that below this value, the total number of topplings scales as O(nlogn), 
whereas for p above this value the total number of topplings increases expo- 
nentially in n. 

Example 7. Consider the fixed-energy Manna model on the complete graph 
K n , with self- loops. A site is active if there are at least two particles, and 
in a toppling each of the two particles chooses a new position uniformly at 
random from its neighbors, which in this case means from all sites. 

We can view this sandpile model as an Ehrenfest urn model, with a 
stopping criterion. We say that the n sites of K n represent n balls in an 
urn, of two colors. One color (green) represents a site with an odd number 
of particles, and one color (red) a site with an even number of particles. In 
an Ehrenfest urn model, a move consists of picking one ball at random from 
the urn, removing it, and instead adding a ball of the opposite color. 

Now consider a toppling of an active site in the sandpile model. We 
imagine that we first (temporarily) remove the two particles from the active 
site, and then one by one add them to their new positions. Removing the 
particles does not change any of the parities. Putting back one particle 
changes the parity of one site chosen uniformly at random from the graph, 
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and likewise for the other particle. Thus, a toppling of an active site is the 
equivalent of two moves in the urn model. 

The sandpile model is stable if and only if each particle is on a different 
site. Therefore, if the particle number is m, then the stopping time for the 
according urn model is the first time that the number of green balls is m (if 
m > n, then the model never stabilizes). 

Now suppose that n is large. Then typically, the number of particles is 
close to np. In the initial configuration, typically the number of sites with an 
odd number of particles, is close to |(1 — e~ 2p ) (n times the probability that 
a Poisson(p) number is odd), which is strictly less than np. In an Ehrenfest 
urn model, from any starting point the fraction of green balls typically tends 
to its equilibrium value of 1/2. 

Therefore, if p < 1/2, then the model meets the stopping criterion while 
tending to equilibrium. It is known that the expected total number of moves 
to reach equilibrium for the first time, starting from no green balls, is of the 
order n log n [22] , only slightly more than linearly increasing with n. However, 
when p > 1/2, then the model does not meet the stopping criterion while 
tending to equilibrium. On the contrary, it meets the stopping criterion only 
in a large deviation from equilibrium. This will typically take an exponen- 
tially large number of moves. Meanwhile, the number of green balls fluctuates 
round its equilibrium number, which we can describe as quasi-stationary be- 
havior. Therefore, there is a sharp transition in the total number of topplings 
needed to stabilize, at p = 1/2. If we monitor the quasi-stationary behavior, 
then we find 1/2 as an estimate for p tr . Actually, the threshold density is 
Pth = 1. 

To find the stationary density, we say that one site is the sink. Then both 
topplings and additions are equivalent to one or more moves in the Ehrenfest 
urn model, this time without stopping criterion. Therefore, p s = 1/2. We 
conclude that this model behaves as conjectured for the Manna model on Z. 

5 Proofs 

5.1 Proof of Theorem Q] 

Our main tool for the proof of Theorem [T] the following lemma. Consider the 
set A = {rj : t](x) = 1 for all < x < n}, which contains two elements which 
we denote by rjf — (1, 1, . . . , 1) and r) e — (1, 1, . . . , 1, 0). 

Lemma 8. Consider the driven sandpile model on I n with toppling environ- 
ment ( given by ((0) = r and ((x) = s for all other x G I n . Then the set 
A defined above is absorbing. More precisely, denoting the addition site by 
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a, transitions between rjf and rj e are as follows. If n — a is even, then the 
configuration flips from r)f to rj e or vice versa; if n — a is odd, then the con- 
figuration does not change upon addition at a and subsequent stabilization. 

Since obviously, A is reachable from any configuration, the lemma implies 
that there are only two recurrent configurations: the 'full configuration' r)f 
in which all sites have height 1, or the 'non-full configuration' t] e in which all 
sites have height 1, except site n. 

Proof. We use the well known description [29] which gives the result of one 
addition and subsequent stabilization in a particle configuration rj in the 
driven model on I n : Let a be the addition site. Let % be max{— l,i < a : 
r)(i) = 0}, that is, % is the closest empty site to the left of a, and if there is 
no such site, then i = —1. Similarly, let j be min{n + l,j > a : = 0}. 
The resulting configuration after the addition and subsequent stabilization 
is denoted £. Then only sites in (i, j) have toppled, and £ is as follows: sites 
i and j have gained 1 particle and site i+j — a has lost 1 particle. Note that 
if 77(a) = 0, then the above description simply states that site a will become 
full. If 77(a) = 1, then at most three sites change height. 

Now we return to the sandpile model with £(0) = r, and let 77 G A 
We will show that if n — a is odd, then we end up in 77 after addition and 
subsequent stabilization. If n — a is even however, then we flip to the other 
element in A. 

We first note that if the particle configuration is non-full and a = n, 
then it is obvious that the result is the full configuration. Now suppose that 
a is a full site. After the addition, we stabilize in steps. In the first step, 
we perform only s-topplings (that is, we do not topple site 0) until all sites 
except site are stable. This step is equivalent to stabilization in the driven 
classical sandpile on [1, n], with sinks at sites and n + 1. In the second step, 
we topple site once, and subsequently again only perform s-topplings until 
all sites except site are stable. We repeat the second step until the whole 
configuration is stable. 

To find the particle configuration after the first step, we apply the above 
description with i — and j = j\ where ji — n if 77 = T) e , and j\ — n + 1 if 
77 = rjf. Then after the first step, site now has two particles, and all other 
sites are full except for site ji — a, which is now empty. 

To find the particle configuration after the second step, we apply the 
above description twice, because toppling site is equivalent to emptying 
site and making two additions to site 1. Therefore, after the second step 
site has again two particles, and all other sites are full except site j\ — a — 2, 
which is now empty. 
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We repeat the second step k times, with k such that either j\ — a — 2k = 
or ji — a — 2k = 1. In the first case, we have reached a stable configuration, 
namely rjf. This will be the case if ji = n and n — a is even, or ji — n + 1 
and n — a is odd. In the second case, all sites will topple one more time, so 
that the final configuration is r] e . This will be the case if ji = n + 1 and n — a 
is even, or j\ = n and n — a is odd. This finishes the proof. □ 

The idea of the proof of Theorem [T] in words is as follows: we think of 
almost the entire state of being built up of 'blocks' in which the leftmost site 
is r, and all other sites are s. Each of these blocks is like the model in Lemma 
|8j after a transient period, either all sites are full, or only the rightmost site 
is empty. Moreover, we will find that at stationarity, these configurations are 
equally likely. Therefore, since the proportion of r sites is f r , the proportion 
of empty sites at stationarity is f r /2. 

Proof of Theorem^ Let (r/, () be a state on I n of the quenched directed 
model. Let R n be the number of r sites in ( n , and call X i} i — 1, . . . , R n the 
positions of the r sites, in increasing order. We define B = [0,X\ — 1] and 
E>i = [Xi, Xi + \ — 1] for i — 1, . . . , R n . We refer to the B^s as a 'block' and 
denote by U be the number of sites in Bi. 

The block B behaves as a classical sandpile, implying that the only 
recurrent (sub) configurations on B are those with at most one site with 
height 0. For all other blocks, Lemma [8] implies that there are two recurrent 
sub-configurations, namely all U sites are full (a 'full' block), or all sites are 
full except the rightmost site X i+ i — 1 (a 'non-full' block). In fact, according 
to Lemma [8j an addition to the a-th site (counting from the left) in B^ with 
U — a odd will leave the configuration unchanged, whereas an addition to 
the a-th site in Bi with Zj — a even will flip the block, that is, change its 
configuration from full to non-full or vice versa. 

After a random but a.s. finite time, all blocks Bi, ... , Br u are either full 
or non-full, and in Bo there is at most one empty site. We will now argue 
that for each block Bi, . . . , -Br„, the probability that it is full converges to 
1 /2 as time tends to infinity. 

Consider the effect of one addition. Suppose we add at a site in Bi for 
some % > 1. As noted before, either Bi flips, or it does not, depending only 
on the location of the addition and its relative position in Bi. In particular, 
whether or not the block flips does not depend on the configuration. If Bi 
flips, then either zero or two particles leave Bi. They have the effect of an 
even number of additions to site Xi + \. Since one addition at Xi + i has either 
no effect, or results in a flip, two additions have the effect of no change at all. 
Therefore, the configuration of block -Bj+i does not change, and as a result, 
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two or zero particles leave B i+1 into B i+2 . Continuing this way, we conclude 
that if Bi flips, then no other block flips. If Bi does not flip, then one particle 
leaves Bi, and it has the effect of one addition to site X i+1 . Therefore, B i+ i 
will flip if U is odd, otherwise, one particle leaves -Bj+i and has the effect of 
one addition to site X i+2 . We conclude that if Bi does not flip, then 'the first 
odd block downstream' flips. In other words, block Bj flips, where lj is odd, 
and there is no i < k < j such that 4 is odd. If there exists no such j, then 
no block flips, so that the addition did not change the configuration. 

In the interval [0,Xi — 1] there is at most one empty site, but it is not 
necessarily site Xi — 1. If the addition site is in B , then the new configuration 
depends on the current configuration. However, at most two particles leave 
this interval: either one site flips from to 1, one site flips from 1 to 0, or 
two sites flip (one from to 1 and one vice versa). In the first two cases, 
zero or two particles leave B Q , so that the configuration downstream is not 
affected. In the third case, one particle leaves B , so that the first odd block 
downstream flips. 

Thus we can partition all sites, except the sites in [0,Xi — 1], into deter- 
ministic non-empty sets Si, i — 1, ... , R n , such that if there is an addition 
to a site in Si, then Bi and only Bi flips. We can summarize the effect of 
an addition as follows: pick a random block, such that for alH = 1, . . . , R n , 
block Bi is chosen with probability j^j-, and flip it. With probability 
flip either a site in [0, Xi — I], or two in [0,Xl — I] and the first odd block. 
With probability ( " +1 ~^ 1 + ~ EJS ' l) , flip nothing. 

It is clear from the above description that for each block B\, B 2 , . . . , -Br„, 
the probability that it has flipped an even number of times until time t 
(inclusive, say) converges to 1/2 as t — > oo. Hence, the expected number 
of sites with height in block B i: converges to 1/2 as t — >■ oo, for all % = 
1, . . . , R n . Block B can be viewed as a classical sandpile model on Xi sites 
(with a slower addition rate of course). Therefore, the expected number of 
sites in B with height converges to l/(x + 1) as t — > oo. 

Writing P for // x v and E for the corresponding expectation, we now 
write 

Pn = ^^Y< E \ibv t W\ R n = k,X 1 = x\F(R n = k,X 1 =x), 
k,x \ i=0 / 
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which is a finite sum. When t — > oo, this converges to 



lim p, 



n 




2(n+l) 



A- 



(x + l)(n + l) 



1 



F(R n — k,X 1 — x) 



1 



E(-R„) 




P(i? n = fc, Xi = x). 



1 



2(n + l) 



Since ^-E(i2 n ) tends to z/(£ : £(0) = r) = / r as n — > oo and the last term 
tends to 0, we obtain 



5.2 Proof of Theorem [2] and Theorem [3] 

Proof of Theorem [#| Part (a) follows from Theorem [TJ hence we need only 
to prove (b). Recall that the initial configuration rj consists of a Poisson-p 
number of particles at every site. It will be convenient to also consider the 
parity of the particle numbers at each site. To this end, let rf = rj mod 2, 
and call r)'(x) the parity of x. We define E{x) = (r](x) — r/(x))/2, and call 
E(x) the number of excess pairs of particles at x. The E(x) will play a crucial 
role in what follows. When a site x is unstable, it topples, which means that 
it loses two particles. If ((x) = r, the two particles can be thought of as 
staying together, while when ((x) = s, the two particles are separated, since 
they are sent to different neighbors. It is for this reason that we need to 
analyze in detail how topplings of an s-site take place, since only around 
such s-sites it is possible that the excess pairs can be absorbed. 

Consider a site x with ((x) = s. For the time being, we restrict our 
attention to {x — 1, x, x + 1}. We topple unstable sites in {x — 1, x, x + 1}, 
until we reach the final stable configuration £ on {x — 1, x, x + 1}. Since we 
work in the restricted model, we have that ({x — 1) = ((x + 1) = r. First, 
observe that topplings will never change the parity of x. 

Furthermore, the sum of the parities of x—1 and x+1, computed mod 2, 
does not change either under topplings. Finally, note that if x topples at all, 
then £(x — 1) can not be 0, and hence £(a; — 1) = 1. Hence, we have the 
following table, dealing with the parities before versus parity/configuration 
after stabilization of {x — l,x,x + 1}, under the assumption that x topples 
at least once. 



lim lim p 1 



1 " fr/2. 



□ 
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parity before stabilization 


configuration after stabilization 


000 


101 


001 


100 


010 


111 


100 


100 


Oil 


110 


101 


101 


110 


110 


111 


111 



It follows from the table that the set {x — l,x, x + 1} can absorb an excess 
pair if and only if the original parities are 000 or 010. We call a site x with 
= s a hole, if rf(x — l) = rf(x + l) = 0. In all other cases, all excess pairs 
that are present in, or enter {x — l,x,x + 1} from the left, will after some 
topplings in {x — 1, x, x + 1} result in excess pairs leaving at the right of the 
set. For every excess pair, this requires at most two topplings per site. This 
observation suggests that in order to decide whether or not a configuration 
is stabilizable, we need to compare the spatial density of excess pairs to the 
density of holes. Or, in other words, we need to compare the probability that 
a given site a; is a hole with the expected number of excess pairs at x. Below, 
we formally prove that this is the correct picture, but first we compute these 
quantities. 

It is easy to compute the (conditional) probability for a site i to be a 
hole, given that ((x) = s. Indeed, since the original configuration is product 
measure with Poisson-p marginals, the probability that a site has parity is 
1(1 + e~ 2p ). It follows that the probability of a site x being a hole is 



fr, 



1 + e 



-2p\2 



(5.1) 



A similar and related computation shows that the expected number of excess 
pairs is equal to 

l(p-kl + e- 2p ))- (5-2) 



Comparison of (5.1) and (5.2) then finishes the proof. 

It remains to show that it is indeed the case that the relation between 
the probability of a hole and the expected number of excess pairs determines 



whether or not stabilization occurs. To this end, we first assume that (5.1) 
< (5.2), that is, the spatial density of holes is smaller than the expected 



number of excess pairs per site. Suppose that stabilization occurs under 
this assumption. First we note that if the configuration stabilizes, all excess 
pairs must be absorbed by a hole. Indeed, if there is an excess pair that is 
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never absorbed, with positive probability, then by ergodicity there must be 
infinitely many such pairs to the left of the origin, which implies that the 
origin will topple infinitely many times. However, a pairing in which every 
excess pair is paired to a unique hole is only possible if the spatial density 
of the holes is at least as large as the spatial density of the excess pairs; this 
follows, for instance, from the well known mass transport principle, see e.g. 
[TO] . This violates the assumption, and we are done. 



For the other direction, we assume that (5.1) > (5.2). Define a random 
walk S(x),x G Z as follows. First we put S = 0. Writing h(x) = 1 if there 
is a hole at x, and h(x) = otherwise, we put 

S(x) =S(x-l) + E(x)-h(x). 

By assumption, S has negative drift. We say that x is extremal if S(y) > S(x) 
for all y < x. The collection of extremal points has a well defined spatial 
density. In particular, there are infinitely many extremal points in both 
directions. It is now simple to see that any excess pair will be absorbed before 
reaching the first extremal point to the right. In particular, no excess pair will 
travel across an extremal point. From this it follows that the configuration 
stabilizes. □ 

Proof of Theorem [3| The proof of this result is very similar to the proof of 
Theorem [2} Again, we have to compare the number of excess pairs to the 
number of holes. For the same reason as before, each hole can accommodate 
one excess pair. What follows from this is that p t h has the same numerical 
value as p tr . □ 



5.3 Proof of Theorem |4] 

The main difficulty of the proof is that we work in infinite volume, so we 
cannot use as a starting point that the dynamics becomes periodic. 

The cases p < 1 and p > 2 are not difficult. The case 1 < p < 2 however 
is rather tedious. Our strategy is to generalize the well-known notion of a 
forbidden sub- configuration (FSC) [I31[25], special local configurations which 
do not appear in the stationary state of the system. We will show that in 
infinite volume with parallel toppling order, there are a number of additional 
local configurations which, essentially, do not appear in the stationary state 
of the infinite-volume system in the sense that their spatial density tends to 
as t — > oo. Once we have proved that, it is easy to show that in all remaining 
possibilities, all sites have period 2, and this will conclude the proof. 

In the proof, the following lemma plays an important role. In [S] it was 
stated and proved for Z 2 , but their proof works for general graphs. We state 
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a version suitable for our purposes. We say that two configurations rj and £ 
are mirror images of each other, if rj(x) = 3 — £(x) for all x. 

Lemma 9 (|8J). Let two height configurations rj and £ be mirror images of 
each other. Then after performing, for each, a parallel toppling time step, 
the two resulting configurations are again mirror images of each other. 

Proof of Theorem [7} The case p < 1 follows from the fact that the system is 
(a.s.) stabilizable when p < 1. This is a well-known fact, proved in several 
papers, see for example [21]. We remark that in |25j . it is moreover proved 
that for p = 1 the system is not stabilizable. 

Next, consider the case p > 2. We have that E p (3 — ?y(0)) < 1. Therefore, 
when rj is chosen according to p p , the mirror image fj of rj is a.s. stabilizable. 
So when we start the dynamics from both rj and fj, then for every site x, 
ff{x) G {0,1} eventually. According to Lemma [9j rf{x) G {2,3} eventually, 
proving this case. (Note that when a site x, and its two neighbours, are all 
unstable, then after toppling all unstable sites, the height of x is the same as 
before. So it is possible to have a well defined limiting configuration, despite 
the fact that topplings continue to take place.) 

Suppose now that 1 < p < 2. In this case, when we choose rj according 
to p p , neither rj nor fj is stabilizable a.s. Therefore, again using Lemma |9j 
for every site x, there is a (random) time T < 00 such that if t > T, x has 
been stable for at least one time step, and also has been unstable for at least 
one time step. Hence, for any finite subset W of Z, there is a (random) time 
T w < 00 such that before time T w , every site in W has been both stabel and 
unstable for at least one time step. It is well known (see [H] and also [25] ) 
that this implies that for t > TV, the total number of particles in W is at 
least equal to the number of internal edges in W. A local configuration in W 
which does not satisfy this property is called a forbidden sub- configuration 
(FSC). It also follows from [TJ] that both rj t {x) and fjt(x) are at most 3, for 
all x G W and t > T w . 

Comparing the dynamics from both rj and its mirror fj, we have from 
Lemma [9] that if rfix) is unstable, then ff(x) is stable, and vice versa. Since 
fj eventually has no FSC in W, we see that in rj, also mirrors of FSC's no 
longer appear eventually. For instance, since 00 is a FSC, and therefore does 
not appear anymore eventually, neither does e.g. 33 occur in the long run, 
on any given position, and similarly for all other FSC's. 

Apart from the absence of FSC's, we claim that also the patterns 01, 10 
and 111 are very rare, in the sense that the spatial density of these finite- 
dimensional patterns converges to as t tends to infinity. Using Lemma [9] 
or the above this implies that also 23, 32 and 222 have small spatial density 
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for t large. We prove this claim at the end of this proof, but first show that 
it implies the desired result. 

It is easy to check from the absence of the above mentioned local configu- 
rations, that apart from local perturbations which occur with spatial density 
tending to 0, the only possible configurations are those in which we alternat- 
ingly choose elements from {0, 1, 11} and {2,3,22}. Any finite-dimensional 
configuration which is built up from alternatingly choosing elements from 
{0, 1, 11} and {2, 3, 22} has period 2. If the full configuration would be like 
that (that is, if the FSC's, 01, 10, 111, 222, 23 and 32 would not occur at all) 
than it would follow that 



a p (t) + a p (t + 1) = 1, (5.3) 

for all t. Indeed, if that were the case, then the configuration would be 
periodic with period 2, with each site alternating between a stable and an 
unstable state. However, we have to deal with the rare perturbations of the 
FSC's and the configurations 01, 10, 111, 222, 23 and 32. This however is not 
a real complication. For any e > 0, we find T so large that the sum of the 
spatial densities of all the problematic sub-configurations is at most e, for all 



t > T. It is clear that for these values of t, the sum in (5.3) is close to 1. 
This is enough to prove the result. 

It remains to prove the claim that the patterns 01, 10 and 111 have spatial 
densities which converges to as t — > oo. We start with 01 and 10. We write 
fi P for the probability measure describing the state of the process at time t, 
hence /i° = /i p . Since the initial measure /z° is ergodic, and /i* is a factor 
of /j,° p , it follows that also /x* is ergodic. Hence both 01 and 10 have a well- 
defined spatial density at time t, which is an a.s. constant, depending on t. 
By symmetry, these two densities must be the same. (We do not need that 
they are the same, but we do need that if one is positive, so is the other.) 

Now let, as before, W be a finite interval, and suppose that for some t, the 
configuration in W contains the patterns 01 and 10 (in that order). When 
rf(x) = and rf(x + 1) = 1, we will say that the pattern 01 occurs from x 
(the leftmost vertex of the pattern) and similarly for other local patterns. 

Suppose 01 occurs from x at a certain time s, where s is much larger than 
Tw, how large exactly will be decided later. We are interested in the last 
time before time s that the configuration at (x, x + 1) was different from 01; 
we call this the previous configuration. So the previous configuration is not 
necessarily the configuration at time s — 1. Instead, it is the last configuration 
at (x,x + 1) which was different from the current one. The at a; can only 
have been obtained through a toppling at x, and since x + 1 has only one 
particle, site x + 2 must have been stable. Hence the previous configuration 
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from x must have been either 200 or 201. Since 00 is an FSC, assuming 
that the previous configuration is at a time still larger than T w , 00 is not 
allowed, so 201 is the only option. This means that the pattern 01 cannot be 
created after time Tw'- if it occurs from x, then the previous configuration 
from x must have been 201, implying that 01 occurs from x + 1. A similar 
(symmetric) argument is valid for the occurrence of 10. 

It follows that when we go backwards in time, as long as we remain larger 
than Tw, the pattern 01 'travels' to the right, and the pattern 10 'travels' 
to the left. This implies that at some time in the past they must have met, 
and we must have either the pattern 0110 or the pattern 01 * 10 in W, where 

* is arbitrary. All this, again, is only valid if we assume that when we go 
backwards in time, we never get to a time smaller than T before one of the 
two patterns 0110 or 01 * 10 occurs. It is clear that we can take s so large 
that this is the case with overwhelming probability, since the model does not 
stabilize. 

But now observe that 0110 is an FSC, and therefore this is not a possibil- 
ity, since FSC's do not occur in W after time Tw- Furthermore, the previous 
configuration of 01*10 is an FSC whenever the * is or 1, and therefore these 
are not possible. If * = 2, then the 2 must have arisen from the previous 
configuration through a toppling of (at least) one of its neighbours, but then 
the next to the toppled site is impossible. Similarly for the possibility that 

* = 3. We conclude that all cases lead to impossibilities, and therefore, on 
any finite interval W, we cannot have the patterns 01 and 10 together, after 
sufficiently long time. This implies that the (common) density of 01 and 10 
must tend to zero, as time tends to infinity. 

Next, we focus on the pattern 111, and again we ask the question what the 
previous configuration could have been. (Recall that the previous configura- 
tion is not the configuration one time unit earlier, but the last configuration 
on that position that was different from 111.) If the previous configuration 
contains a 0, then exactly one of the neighbours of this must be unstable. 
Indeed, after the toppling of all unstable sites, this is the only way to go 
from a to a 1. But this implies that the other neighbour of the is a or 
a 1, which means that the previous configuration contains either an FSC, or 
the configuration 01 (or 10), which we know by the previous paragraph has 
density tending to 0. 

If the previous configuration contains a 1, then the neighbours of this 1 
can not be unstable, and by the same argument as above, the neighbours of 
this 1 cannot be or 1 (after time Tw)- Hence, the only possibility is that 
the previous configuration would be 111, but this is a contradiction, since 
by definition, the previous configuration must be different from the current 
one. It remains to rule out all previous configurations which consist of only 
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2's and 3's. It is obvious however that it is not possible to reach 111 through 
such a previous configuration: if all sites are unstable then the central site 
receives two particles and is again unstable. 

We summarize that if the pattern 111 occurs from i at a given time 
larger than Tw, then either it was present from x at time TV or one of 
the patterns 01 or 10 was present in the previous configuration. The first 
possibility requires that the central site has not toppled since time T w , which 
has probability tending to as time tends to infinity. The second possibility 
gives us that the spatial density of the pattern 111 tends to 0, just as we 
showed for the patterns 01 and 10. □ 

Acknowledgment: We thank Lionel Levine, David Wilson and Alexander 
Holroyd for stimulating discussions. 
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